pop_size
## [1] 500
number_of_generations
## [1] 800
gs_vec
## [1] 3 7 13 18 41
fnc_vec
## [1] 1
initial_variance_vec
## [1] 0
initial_knowledge_vec
## [1] 0
average_pdrift_vec
## [1] 0.1 0.2 0.3
variance_pdrift
## [1] 0
time_cost_vec
## [1] 0.01 0.02 0.03
collective_vec
## [1] 0 1
result_grid
#Static environment with evolving sharing and innovation
load(paste0(name,"/",1,name,".RData"))
final_data_all<-expand.grid(Group_Size=gs_vec,FNC=fnc_vec,average_pdrift=average_pdrift_vec,initial_knowledge=initial_knowledge_vec,initial_variance=initial_variance_vec,time_cost=time_cost_vec,collective=collective_vec,m1=NA,sd1=NA,m2=NA,sd2=NA,m3=NA,sd3=NA,m4=NA,sd4=NA,m5=NA,sd5=NA,m6=NA,sd6=NA,m7=NA,sd7=NA)
final_data_all[1,]<-final_data[1,]
for (loop_parameters in 1:nrow(final_data)){
if( file.exists(paste0(name,"/",loop_parameters,name,".RData")) ){
load(paste0(name,"/",loop_parameters,name,".RData"))
final_data_all[loop_parameters,]<-final_data[loop_parameters,]
} else {
print("File missing")
}
if (runif(1)<1.2){
par(mfrow=c(2,4))
matplot((result[[1]]),type="l",xlab="Generation",ylab="bound_initial_bias",col="gray",ylim=c(lb[1],ub[1]));lines(1:number_of_generations,apply(result[[1]],1,mean),lwd=2)
matplot((result[[2]]),type="l",xlab="Generation",ylab="bound_theta",col="gray",ylim=c(lb[2],ub[2]));lines(1:number_of_generations,apply(result[[2]],1,mean),lwd=2)
matplot((result[[3]]),type="l",xlab="Generation",ylab="delta",col="gray",ylim=c(lb[3],ub[3]));lines(1:number_of_generations,apply(result[[3]],1,mean),lwd=2)
matplot((result[[4]]),type="l",xlab="Generation",ylab="SI use",col="gray",ylim=c(lb[4],ub[4]));lines(1:number_of_generations,apply(result[[4]],1,mean),lwd=2)
matplot((result[[5]]),type="l",xlab="Generation",ylab="First",col="gray",ylim=c(lb[5],ub[5]));lines(1:number_of_generations,apply(result[[5]],1,mean),lwd=2)
matplot((result[[6]]),type="l",xlab="Generation",ylab="Perfomance",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[6]],1,mean),lwd=2)
matplot((result[[7]]),type="l",xlab="Generation",ylab="Perfomance of first",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[7]],1,mean),lwd=2)
Sys.sleep (.2)
}
}
#Static environment with evolving sharing and innovation
load(paste0(name,"/",1,name,".RData"))
final_data_all<-expand.grid(Group_Size=gs_vec,FNC=fnc_vec,average_pdrift=average_pdrift_vec,initial_knowledge=initial_knowledge_vec,initial_variance=initial_variance_vec,time_cost=time_cost_vec,collective=collective_vec,m1=NA,sd1=NA,m2=NA,sd2=NA,m3=NA,sd3=NA,m4=NA,sd4=NA,m5=NA,sd5=NA,m6=NA,sd6=NA,m7=NA,sd7=NA)
final_data_all[1,]<-final_data[1,]
for (loop_parameters in 1:nrow(final_data)){
if( file.exists(paste0(name,"/",loop_parameters,name,".RData")) ){
load(paste0(name,"/",loop_parameters,name,".RData"))
final_data_all[loop_parameters,]<-final_data[loop_parameters,]
} else {
print("File missing")
}
if (runif(1)<1.2){
par(mfrow=c(2,4))
matplot((result[[1]]),type="l",xlab="Generation",ylab="bound_initial_bias",col="gray",ylim=c(lb[1],ub[1]));lines(1:number_of_generations,apply(result[[1]],1,mean),lwd=2)
matplot((result[[2]]),type="l",xlab="Generation",ylab="bound_theta",col="gray",ylim=c(lb[2],ub[2]));lines(1:number_of_generations,apply(result[[2]],1,mean),lwd=2)
matplot((result[[3]]),type="l",xlab="Generation",ylab="delta",col="gray",ylim=c(lb[3],ub[3]));lines(1:number_of_generations,apply(result[[3]],1,mean),lwd=2)
matplot((result[[4]]),type="l",xlab="Generation",ylab="SI use",col="gray",ylim=c(lb[4],ub[4]));lines(1:number_of_generations,apply(result[[4]],1,mean),lwd=2)
matplot((result[[5]]),type="l",xlab="Generation",ylab="First",col="gray",ylim=c(lb[5],ub[5]));lines(1:number_of_generations,apply(result[[5]],1,mean),lwd=2)
matplot((result[[6]]),type="l",xlab="Generation",ylab="Perfomance",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[6]],1,mean),lwd=2)
matplot((result[[7]]),type="l",xlab="Generation",ylab="Perfomance of first",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[7]],1,mean),lwd=2)
Sys.sleep (.2)
}
}
print(simulation)
## function (plot_yes, n, initial_variance, initial_knowledge, average_pdrift,
## variance_pdrift, initial_bias, theta, delta, false_costs,
## time_cost)
## {
## p_drift = rnorm(n, average_pdrift, variance_pdrift)
## base_line = 0.5
## d_time = rep(NA, n)
## h = 0.05
## sigma = 1
## t = 1
## predator = runif(1) < base_line
## if (predator)
## intercept = rnorm(n, 0.5 * theta + initial_knowledge +
## initial_bias, initial_variance)
## if (!predator)
## intercept = rnorm(n, 0.5 * theta + initial_knowledge -
## initial_bias, initial_variance)
## L = intercept
## attribute = 1
## dec = rep(FALSE, n)
## social_strength = 0
## new_dec = TRUE
## if (plot_yes == 1) {
## L_large = L
## }
## while (any(dec == FALSE)) {
## dec_before <- dec
## d_time[!(L <= 0 | L >= theta) == dec] = t
## dec = L <= 0 | L >= theta
## if (new_dec == TRUE) {
## if (sum(dec) > 0) {
## social_strength = (sum(L[dec] >= theta[dec] *
## 0.5) - sum(L[dec] <= theta[dec] * 0.5))
## }
## }
## randomNorm = rnorm(n, 0, 1)
## evidence = (p_drift + (social_strength * delta)) * h +
## (sigma * h^0.5) * randomNorm
## evidence[dec] = 0
## L = L + evidence
## t = t + 1
## new_dec <- any((dec) != (dec_before))
## if (plot_yes == 1) {
## L_large <- cbind(L_large, L)
## }
## }
## d_time[is.na(d_time)] <- t
## d_time = (d_time - 1) * h
## final_decision <- L >= theta * 0.5
## if (predator)
## fitness = ifelse(final_decision, 1, -false_costs[2]) -
## (d_time * time_cost)
## if (!predator)
## fitness = ifelse(final_decision, 1, -false_costs[1]) -
## (d_time * time_cost)
## sim_data <- data.frame(final_decision, d_time, fitness, predator)
## if (plot_yes == 1) {
## matplot(seq(h, t * h, h), t(L_large), type = "l", ylim = c(0,
## max(theta)))
## }
## return(sim_data)
## }
print(runEA)
## function (numAgents, pop_size, number_of_generations, gif, false_costs,
## time_cost, initial_knowledge, initial_variance, average_pdrift,
## variance_pdrift, collective)
## {
## if (collective == 0)
## n_eval = ceiling(pop_size/as.numeric(numAgents) * 25)
## if (collective == 1)
## n_eval = ceiling(pop_size * 10)
## library(scales)
## result = NULL
## init = "dev"
## genpool = matrix(NA, pop_size, length(ub))
## if (init == "rand") {
## for (gene_i in 1:length(ub)) {
## start_mean_gen <- runif(1, 0.4, 0.6)
## a = 1
## b = (a/start_mean_gen) - a
## chromosom = rbeta(pop_size, a, b)
## chromosom = rescale(chromosom, c(lb[gene_i], ub[gene_i]),
## c(0, 1))
## genpool[, gene_i] <- chromosom
## }
## }
## else {
## inits <- c(0, 1, 0.1, 0, 0)
## for (gene_i in 1:length(ub)) {
## genpool[, gene_i] <- inits[gene_i]
## }
## }
## result = rbind(result, c(apply(genpool, 2, mean), 0, 0))
## fitness_long = NULL
## results_long = NULL
## for (gen in 1:(number_of_generations - 1)) {
## probs = rep(1, pop_size)
## fitness = rep(0, pop_size)
## count = rep(0, pop_size)
## first_perf = rep(NA, n_eval)
## for (evals in 1:n_eval) {
## if (collective == 0) {
## index = sample(1:pop_size, numAgents, prob = probs)
## }
## else {
## index = sample(1:pop_size, 1, prob = probs)
## }
## count[index] = count[index] + 1
## probs[index] = probs[index]/5
## probs <- probs/max(probs)
## plot_yes = 0
## deltas = genpool[index, 3]
## deltas[genpool[index, 4] == 0] = 0
## deltas[genpool[index, 5] == 1] = 10
## if (collective == 0) {
## sim_result <- simulation(plot_yes, numAgents,
## initial_variance, initial_knowledge, average_pdrift,
## variance_pdrift, initial_bias = genpool[index,
## 1], theta = genpool[index, 2], delta = deltas,
## false_costs, time_cost)
## fitness[index] <- fitness[index] + sim_result$fitness
## }
## else {
## sim_result <- simulation(plot_yes, numAgents,
## initial_variance, initial_knowledge, average_pdrift,
## variance_pdrift, initial_bias = rep(genpool[index,
## 1], numAgents), theta = rep(genpool[index,
## 2], numAgents), delta = rep(deltas, numAgents),
## false_costs, time_cost)
## fitness[index] <- fitness[index] + mean(sim_result$fitness)
## }
## firsts = sim_result$fitness[sim_result$d_time ==
## min(sim_result$d_time)]
## means = mean(firsts)
## first_perf[evals] <- firsts[nnet::which.is.max(-abs(firsts -
## means))]
## }
## if (gif == 1) {
## fitness_long = c(fitness_long, (fitness))
## results_long = rbind(results_long, (genpool))
## }
## k = 5
## winners1 = winners2 = rep(NA, pop_size)
## for (turns in 1:pop_size) {
## indexes = sample(1:pop_size, k)
## winners1[turns] = indexes[which.max(fitness[indexes])]
## }
## genpool <- genpool[winners1, ]
## for (cross in 1:ncol(genpool)) {
## corssover_rate = 0.01
## dummy_index <- which(runif(pop_size) < corssover_rate)
## dummy_index_partner <- sample(pop_size, length(dummy_index))
## if (length(dummy_index) > 0) {
## for (corss_process in 1:length(dummy_index)) {
## which_gen <- sample(length(ub), 1)
## if (which_gen %in% c(4, 5))
## which_gen <- c(4, 5)
## genpool[c(dummy_index[corss_process], dummy_index_partner[corss_process]),
## which_gen] = genpool[c(dummy_index_partner[corss_process],
## dummy_index[corss_process]), which_gen]
## }
## }
## }
## for (mut in 1:ncol(genpool)) {
## mutation_rate = 0.01
## dummy_index <- runif(pop_size) < mutation_rate
## if (mut <= 3) {
## genpool[dummy_index, mut] <- genpool[dummy_index,
## mut] + rnorm(sum(dummy_index), 0, (ub - lb)[mut] *
## 0.05)
## genpool[genpool[, mut] < lb[mut], mut] <- lb[mut] -
## (genpool[genpool[, mut] < lb[mut], mut] - lb[mut])
## genpool[genpool[, mut] > ub[mut], mut] <- ub[mut] -
## (genpool[genpool[, mut] > ub[mut], mut] - ub[mut])
## }
## else {
## genpool[dummy_index, mut] = rbinom(sum(dummy_index),
## 1, 0.5)
## if (mut == 4)
## genpool[dummy_index, mut + 1][genpool[dummy_index,
## mut] == 1] = 0
## if (mut == 5)
## genpool[dummy_index, mut - 1][genpool[dummy_index,
## mut] == 1] = 0
## }
## }
## result = rbind(result, c(apply(genpool, 2, mean), ifelse(collective ==
## 0, mean(fitness)/(n_eval/(pop_size/numAgents)), mean(fitness)/(n_eval/(pop_size))),
## mean(first_perf)))
## }
## if (gif == 1) {
## library(gganimate)
## library(tweenr)
## library(magick)
## library(purrr)
## library(tidyr)
## library(dplyr)
## library(cowplot)
## plot_data = data.frame(ease = rep("cubic-in-out", length(fitness_long)),
## fitness_long, x1 = results_long[, 1], x2 = results_long[,
## 2], x3 = results_long[, 3], x4 = results_long[,
## 4], x5 = results_long[, 5], id = rep(1:pop_size,
## (number_of_generations - 1)), gen = sort(rep(1:(number_of_generations -
## 1), pop_size)))
## number_of_frames = number_of_generations * 2
## if (number_of_generations > 99) {
## number_of_frames = round(number_of_generations/2)
## plot_data = plot_data %>% filter(gen %in% seq(1,
## number_of_generations, 5))
## }
## plot_data2 = plot_data %>% gather(type, phen, c(x1, x2,
## x3, x4, x5))
## plot_data2$max = ifelse(plot_data2$type == "x1", ub[1],
## ifelse(plot_data2$type == "x2", ub[2], ub[3]))
## plot_data2$min = ifelse(plot_data2$type == "x1", lb[1],
## ifelse(plot_data2$type == "x2", lb[2], lb[3]))
## plot_data2 = plot_data2 %>% mutate(type = ifelse(type ==
## "x1", "Bias", ifelse(type == "x2", "Boundary Seperation",
## ifelse(type == "x3", "Social Drift", ifelse(type ==
## "x4", "SI use", "Follow")))))
## pnew = ggplot(plot_data2, aes(x = phen, y = fitness_long,
## color = as.factor(id))) + geom_point() + transition_states(gen,
## 1, 1) + facet_wrap(~type, scales = "free") + theme(legend.position = "none") +
## labs(title = "Generation: {closest_state}", x = "Phenotype",
## y = "Fitness") + ylim(0, max(plot_data2$fitness_long))
## anim_save(paste("Output/tournement", "test", ".gif"),
## animation = animate(pnew, nframes = number_of_frames,
## width = 1100, height = 600))
## }
## cbind(result)
## }